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Rayleigh-Taylor  Instability  in  Compressible  Media 


David  L.  Book 

Laboratory  for  Computational  Physics 
Naval  Research  Laboratory 
Washington,  D.C.  20375 


Introduction 

The  Rayleigh-Taylor  instability 1-2  occurs  when  a  fluid  supports  a  denser 
fluid  against  gravity,  whereupon  the  two  tend  to  interchange  positions.  It 
is  encountered  frequently  in  nature  and  in  the  laboratory.  For  example, 
inertial  confinement  fusion  experiments,  in  which  an  ablatively  driven  medium 
implodes,  compressing  the  material  ahead  of  the  ablation  front  to  high 
densities,  can  exhibit  Rayleigh-Taylor  instabilities  in  the  ablative  region, 
at  the  compression  front,  or  (in  the  case  of  a  layered  target)  at  an 
interface  between  layers  of  different  density. 

When  the  time  scale  associated  with  the  growth  of  the  instability  is 

short  compared  with  the  time  (kc^)”'1  for  sound  to  traverse  a  wavelength  2r/k, 

one  should  expect  to  have  to  include  the  finite  compressibility  of  the  fluid 

in  calculating  instability  growth  rates.  It  is  not  obvious  a  priori  whether 

finite  compressibility  acts  to  increase  or  decrease  growth  rates.  For 

example,  compression  absorbs  some  energy  that  might  otherwise  go  into  fluid 

motion.  On  the  other  hand,  a  compressible  system  exhibits  more  modes  of 

propagation  than  an  incompressible  one,  so  the  most  unstable  one  might 
Manuscript  approved  February  22,  1984. 
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possibly  have  a  faster  growth  rate  then  the  most  unstable  mode  in  an 
incompressible  medium.  Moreover,  Bernstein  et  al. 8  have  shown  that  in  a 
broad  class  of  general  compressible  hydromagnetic  systems,  the  unstable  modes 
with  lowest  threshold  are  associated  with  incompressible  perturbations.  It 
is  thus  conceivable  that  compressible  and  incompressible  systems  might 
display  the  same  Rayleigh-Taylor  growth  rate. 

Theoretical  research  into  the  Rayleigh-Taylor  instability  can  be  divided 
into  analytical  and  computational  approaches.  Some  of  the  early  analytical 
work  done  on  Rayleigh-Taylor  instability  in  compressible  fluids  was 
inconclusive  or  erroneous.  Vandervoort ^  and  Plesset  and  Hsieh8  both  analyzed 
the  instability  at  the  interface  between  two  polytropic  media.  Recently 
Shivamoggi 6  pointed  out  that  these  two  papers  disagree  with  one  another. 
Replying  to  his  comment,  Plesset  and  Prosperetti7  attributed  the 
contradiction  to  an  error  in  Vandervoort * s  analysis  which  invalidates  the 
latter1 s  treatment  except  when  y  =  1.  They  went  on  to  derive  in  a  simple 
manner  a  general  dispersion  relation  for  an  arbitrary  equation  of  state. 

This  derivation,  however,  itself  makes  use  of  an  identity  which  is  strictly 
true  only  for  y=  1,  namely  the  statement  that  the  ratio  of  the  perturbed 
pressure  to  the  perturbed  density  is  equal  to  the  square  of  the  unperturbed  • 
sound  speed. 

Blake8  wrote  down  without  derivation  a  dispersion  relation  for  Rayleigh- 
Taylor  instability  in  compressible  fluids  which  in  fact  is  correct  for  the 
interface  between  two  isothermal  (uniform-temperature)  fluids  satisfying  an 
isothermal  ( y  =  1)  equation  of  state,  and  argued  that  compressibility  effects 
are  negligible  except  in  the  long-wavelength  limit. 
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Matthews  and  Bluraenthal 3  derived  the  same  "isothermal-isothermal"  formula 

2 

with  the  inclusion  of  volume  radiation  forces  (  a  p  )  .  [The  formally 

identical  dispersion  relation  for  waves  propagating  in  a  medium  consisting  of 

two  stably  stratified  isothermal  layers  is  well  known  to  atmospheric 

scientists;  e.g. ,  Tolstoy  derives  it  in  his  review  article1^  for  the  case  in 

which  the  fluids  have  an  (identical)  arbitrary  y  and  analyzes  the  different 

waves  which  arise.]  McCrory  et  al.11  made  the  physically  plausible  argument 

that  pressure  differences  cannot  be  transmitted  across  the  mode  structure  on 

a  time  scale  shorter  than  (kc^)  \  so  that  compressibility  effects  must  limit 

growth  rates  to  values  <  kcg.  Takabe  and  Mima12  wrote  down  a  variational 

form  for  the  growth  rate  in  the  presence  of  both  compressibility  and  thermal 

conduction,  but  neglected  to  state  clearly  the  assumptions  they  employed. 

Scannapieco 1 3  investigated  the  stability  of  a  slab  of  ideal  polytropic  gas 

with  an  exponentially  increasing  or  decreasing  density  profile  confined 

between  two  rigid  horizontal  walls  separated  by  a  distance  d  and  found  that 

growth  rates  are  enhanced  by  compressibility.  However,  in  analyzing  the 

limit  in  which  the  scale  height  H  satisfies  H  »  d,  he  allowed  the  density  to 

vary  while  treating  the  sound  speed  as  a  constant,  which  is  inconsistent 

unless  the  density  decreases  in  the  upward  direction.  Baker14  solved  the 

Blake 8  dispersion  relation  numerically  for  the  ratio  of  the  square  of  the 

calculated  growth  rate  to  the  incompressible  value  as  a  function  of  Atwood 

2  2 

number  A  =  (  -  p^)/(  0^  +  P^)>  the  ratio  c^/c^  ° ^  t*ie  s<^uares  the  sound 

2 

speeds,  and  the  nondimensionalized  wavelength  g/2kc^.  Because  of  the 
assumptions  implicit  in  the  "isothermal-isothermal"  model,  however,  the  first 
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two  parameters  are  not  independent:  c^/c^  "  ^2^1  ”  O  +  A) /Cl  “  A).  His 

conclusion  that  finite  compressibility  is  sometimes  stabilizing  and  sometimes 

destabilizing  is  therefore  called  into  question. 

On  the  whole,  the  problem  with  the  analytical  approach  to  studying 

compressibility  effects  on  the  Rayleigh-Taylor  instability  has  most 

frequently  been  a  failure  to  specify  clearly  the  model  being  investigated. 

Many  authors  have  attempted  to  derive  model-independent  dispersion  relations, 

or  at  least  formulas  of  wide  applicability,  which  would  not  be  restricted  to 

a  particular  type  of  density  profile.  For  a  medium  to  behave  compressibly 

with  respect  to  a  mode  of  wavenumber  k,  however,  the  dimensionless  wavenumber 
2 

must  satisfy  kcg/g  =  kH  <  1,  where  H  is  the  equilibrium  scale  height.  But  h, 
the  vertical  scale  of  the  perturbation,  typically  satisfies  h  ~  k  \  so  we 
must  have  h  >  H.  The  mode  samples  the  vertical  variation  of  the  equilibrium 
state  and  must  therefore  depend  sensitively  on  it. 

It  is  actually  quite  easy  to  show  for  ideal  polytropic  media  that 
compressibility  destabilizes  the  Rayleigh-Taylor  instability.15  The  energy 
principle  of  Bernstein  et  al. 3  (an  extension  of  the  version  given  earlier  by 
Chandrasekhar16)  predicts  that  polytropic  media  with  finite  adiabatic  index  y 
exhibit  maximum  growth  rates  which  decrease  with  increasing  y. 

Incompressible  fluids,  which  are  obtained  in  the  limit  y  -►  »,  are  thus  more 
stable  then  compressible  ones.  This  has  been  confirmed  experimentally  by 
Asay  (see  Ref.  14). 

Most  computational  approaches  to  the  problem  strive  so  hard  for  realism 
that  they  treat  too  many  physical  processes  simultaneously.  When  something 
happens  in  a  calculation  it  is  hard  to  say  which  process  is  responsible, 
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especially  when  all  the  parameters  in  the  code  have  been  chosen  so  as  to 
simulate  a  particular  laboratory  experiment*  While  such  simulations  are  a 
major  reason  for  computation,  they  are  worthless  unless  every  effect  included 
in  the  model  has  been  carefully  validated  against  analytical  theory  or 
reliable  measurements.  Failure  to  do  this  properly  vitiated  some  early 
simulations  of  Rayleigh-Taylor  instability  in  imploding  pellets.  Another 
difficulty  arose  from  the  necessity  of  performing  well-resolved 
multidimensional  calculations  in  place  of  the  one-dimensional  ones  used  for 
studying  unperturbed  implosions. 

A  way  around  this  difficulty  was  found  using  so-called  ‘’piggyback” 
codes.  17  The  linearized  equations  of  motion  are  analyzed  into  a 
superposition  of  angular  (e.g.,  spherical)  harmonics,  and  the  equations  for 
the  radius-dependent  amplitudes  corresponding  to  one  or  more  such  modes  are 
advanced  in  time  together  with  the  zeroth-order  quantities.  This  technique, 
also  successfully  employed  in  connection  with  incompressible  cylindrical 
liner  implosions,1®  eliminates  the  numerical  resolution  problem.  The  work  of 
Shiau  et  al17  clearly  showed  for  the  first  time  that  flow  of  plasma  across  an 
interface  (a  process  which  can  only  occur  in  compressible  media) ,  while 
stabilizing,  is  not  by  itself  able  to  completely  suppress  the  Rayleigh-Taylor 
instability.  Subsequently,  improvements  in  computational  methods  and 
techniques  of  code  validation  have  enabled  fully  nonlinear  two-dimensional 
calculations  to  be  carried  out  which  predict  the  linear  and  nonlinear 
evolution  of  the  Rayleigh-Taylor  instability  with  high  accuracy. 

Stimulated  partly  by  experiment  and  partly  by  code  results,  quite 
comprehensive  theories  have  now  been  developed  which  take  into  account  such 
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diverse  effects  as  vortex  shedding,  compressibility ,  thermal  conduction,  and 
ablation. 21  In  the  remainder  of  this  paper  I  will  not  be  saying  anything 
further  about  the  use  of  numerical  simulations  to  study  the  Rayleigh-Taylor 
instability. 

A  different  question,  which  I  not  will  also  not  be  considering  here, 
involves  the  time  required  to  establish  a  state  as  a  result  of  an  initial 
localized  disturbance.  This  process,  which  is  instantaneous  in  an 
incompressible  fluid,  lasts  a  time  equal  to  that  required  for  a  sound  wave  to 
propagate  a  few  times  back  and  forth  across  the  entire  system.  Instead,  I 
will  consider  linear  eigenmodes,  which  by  definition  are  initiated  in 
"prepared”  states  involving  the  entire  system.  Finding  the  eigenmodes  in 
compressible  fluids  presents  enough  analytical  difficulty  to  dissuade  one 
from  seeking  the  solution  of  the  general  initial  value  problem.  Relatively 
few  papers  have  been  written  on  this  topic. 

This  paper  is  organized  as  follows.  A  simple  version  of  an  argument 
originally  employed  by  Schwarzschild  (see,  e.g.,  Landau  and  Lifshitz22)  in 
discussing  hydrodynamic  interchange  is  used  to  derive  threshold  criteria  for 
the  Rayleigh-Taylor  and  convective  instabilities  in  arbitrary  stratified 
media.  Then  the  energy  principle  is  used  to  show  that  compressibility  is 
always  destabilizing,  and  Newcomb’s  extension  of  this  result  to  higher 
eigenmodes  of  the  system  is  presented.22  Next,  the  exact  dispersion  relation 
for  the  Rayleigh-Taylor  instability  at  the  interface  between  two  ideal 
polytropic  fluids  with  different  adiabatic  indices,  each  fluid  having  uniform 
temperature,  is  derived  following  Bernstein  and  Book,  214  and  various  limiting 
cases  of  this  result  are  discussed.  The  extension  to  an  arbitary  piecewise 
isothermal  equilibrium  is  sketched,  concluding  the  portion  of  the  paper 

dealing  with  stability  of  static  equilibria. 
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The  only  nonstationary  fluid  states  in  which  the  Rayleigh-Taylor 
instability  can  be  treated  analytically  are  self-similar  expansions  or 
contractions  for  which  the  velocity  is  proportional  to  the  distance  from  the 
center  of  symmetry  (uniform  self-similar  motions).25  Following  a  summary  of 
the  formalism  used  to  discuss  stability  of  such  states,  examples  are  given, 
first  for  implosions  and  then  for  expansions.  A  final  section  summarizes  the 
main  results  and  attempts  to  draw  them  together  by  pointing  out  the  common 
themes  that  run  through  all  of  these  examples. 
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Physical  Basis  of  Gravitational  Interchange  Instabilities 

Suppose  a  fluid  with  vertical  density  and  pressure  profiles  p(y),  p(y) 
is  in  hydrostatic  equilibrium: 


-§7 +  *  ■  °> 


(1) 


where  g  is  the  constant  gravitational  acceleration.  Note  that  p(y)  must 
decrease  monotonically  as  a  function  of  y,  but  p(y)  need  not.  We  assume  that 
the  adiabatic  index  (ratio  of  specific  heats)  y  is  constant. 

Now  consider  an  element  of  fluid  with  differential  volume  AV  at  some 
arbitrary  height  y.  It  contains  mass  An  =  pAV  and  internal  energy  pAV/(y-l). 
Assume  that  it  is  displaced  adiabatically  to  a  new  position  y!,  where  it 
occupies  a  new  volume  AV*  at  a  new  density  p*  and  pressure  p! .  By 
conservation  of  mass, 


PT  AV1  =  pAV  *  Am; 


(2) 


by  adiabatic  invariance  (entropy  conservation), 

p'(  AV')Y=  p(  av)y.  ( 

To  make  room  for  the  displaced  parcel  of  fluid,  a  second  differential 


volume  AV  with  initial  density  p  and  pressure  p  is  displaced  from  location  y 
*=  y*  to  the  first  location  y!  =  y.  We  assume  that  y  -  y  =  h  >  0,  i.e.  ,  the 
first  location  is  above  the  second  (possibly  by  a  finite  distance). 

Evidently  the  second  parcel  of  fluid  after  displacement  has  density  p!  and 
pressure  p*  satisfying 


p’  AV’  =  pAV  =  Am, 

(4) 

"p’  (AV’  )  Y  =  pfAVY). 

(5) 
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Since  the  displacements  are  adiabatic,  the  total  change  in  internal 
energy  is 

flW  =  (pf  AV'  -  pAV  +  A V*  -  pAV)/(Y-l).  (6) 

This  is  accompanied  by  a  total  net  change  in  gravitational  energy  given  by 

6WG  =  An  g(y?  -  y)  +  Am  g(y*  -  y)  =  (  Am  -  Am)gh.  (7) 

If  after  the  interchange  the  displaced  fluids  are  in  equilibrium  with  their 
surroundings  but  the  total  change  in  energy  is  negative  (i.e.,  energy  is 
reduced) , 

6W  =  aw  +  <  0,  (8) 

the  interchange  is  energetically  favored,  and  the  configuration  is  therefore 
unstable . 

There  are  two  conditions  for  a  displaced  differential  volume  to  be  in 
equilibrium  with  its  surroundings:  a  kinematic  condition,  that  it  ’’fill  the 
hole”  left  by  its  counterpart,  and  a  dynamic  conditions,  that  it  be  in 
pressure  balance.  For  the  first  parcel,  these  requirements  imply 

—  (9) 

and 


p'  =  p; 

for  the  second,  they  imply 

AV’  =  AV 


(10) 


(ID 


and 

?'  -  P.  (12) 

Equations  (3)  and  (5)  now  both  reduce  to  a  relation  connecting  AV  and  AV: 
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p  (AV)T=  p  (AV)Y.  (13) 

If  we  calculate  the  energy  in  each  parcel  of  fluid  as  it  undergoes 
displacement  we  see  that  the  work  done  on  the  first  one  by  the  surrounding 
fluid  exactly  equals  the  expansion  work  done  by  the  second.  Thus 

SH  «  0,  (14) 

and  so 

w  =  SWg  =  pAVgh  [p/p)(p/p)  ^  y-l  ] .  (15) 

We  consider  three  cases  of  interest. 


Case  I:  Discontinuous  change  in  density.  Here  we  can  take  y  and  y  to  be  on 
opposite  sides  of  the  discontinuity,  but  contiguous,  so  that  h  is  very 
small.  Since  the  pressure  is  continuous,  AV  =  A V,  and  Eq.  (15)  reduces  to 

W  -  (“p  -  P)  AVgh.  (16) 

The  system  is  thus  unstable  if  p  >  p.  This  is  the  usual  criterion  for  the 
Rayleigh-Taylor  instability. 


Case  II:  Continuous  density  variation,  p  >  p.  Now  y  and  y  need  not  be 

contiguous.  Since  p  >  p  always  holds,  is  again  negative  if  p  >  p.  The 

G 

limit  as  h  ->•  0  yields  the  criterion  for  Rayleigh-Taylor  instability  in  a 
smoothly  stratified  medium, 


Vp»?p  <  0. 


(17) 
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Case  III:  Continuous  density  variation,  p  <  p. 

Even  if  the  density  decreases  monotonically  in  the  upward  direction,  it  is 
still  possible  to  satisfy  <  0,  provided  that 

p/"pT  >  p/  pY.  (18) 

This  is  the  criterion  for  convective  instability. 22 

We  thus  see  that  in  certain  circumstances  a  fluid  stratified  under 
gravitational  acceleration  can  be  unstable  to  overturning,  or  interchange. 
When  the  instability  is  driven  by  a  density  inversion  (dense  fluid  lying 
above  less  dense),  the  name  "Rayleigh-Taylor  instability*'  is  used.  By  means 
of  the  energy  principle,  these  ideas  can  be  carried  a  step  further  to  exhibit 
the  manner  in  which  the  degree  of  compressibility  affects  stability. 
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Energy  Principle 


The  equations  describing  the  evolution  of  a  small  perturbation  about  a 
specific  equilibrium  state  are  linear  in  the  perturbed  fluid  variables.  They 
can  be  reduced  to  a  single  homogeneous  differential  equation  in  terms  of  one 
dependent  variable,  e.g. ,  the  perturbed  pressure  or  displacement.  If  we 
assume  time  dependence  ~  exp  (-iat),  where  a)  is  the  frequency,  an  ordinary 
differential  equation,  usually  of  second  order,  results  for  the  spatial 
dependence.  The  quantity  to  (or  to2)  appears  as  an  eigenvalue,  determined  by 
solving  the  equation  subject  to  appropriate  boundary  conditions.  If  this 
eigenvalue  problem  is  of  Sturm-Liouville  form,  a  number  of  rigorous  theorems 
apply.  The  most  important  of  these  says  that  to  can  be  found  from  a 
variational  principle,  i.e.,  by  looking  for  the  extremal  (usually  minimum) 
value  of  some  functional  over  a  class  of  these  functions  which  satisfy  the 
boundary  conditions  and  other  physical  constraints.  In  hydrodynamic 
stability  problems,  the  variational  principle  has  a  natural  interpretation 
in  terms  of  energies. 

Assuming  an  adiabatic  equation  of  state,  the  potential  energy  W 
associated  with  a  general  small  perturbation  £(x)  about  an  allowed 
equilibrium  of  the  ideal  magnetohydrodynamic  equations  can  be  expressed  in 


the  form- 


w  =  W,+  y  f  d*x  p  (  V*0  2, 


(19) 


where  is  a  quadratic  functional  of  £  which  is  independent  of  y. 
eigenvalues  resulting  from  solution  of  the  Sturm-Liouville  problem 


If  the 


determining  £  are  ordered  by  magnitude  according  to 


'JO2  <  U)2  <  !D2  <  . . . , 


(20) 
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then  the  lowest  (most  unstable)  eigenvalue  is  determined  by  a  variational 
principle 


=  min  (W/K), 


(21) 


where  K  is  a  second  (nonnegative)  quadratic  functional  and  the  minimum  is 
taken  over  all  5  satisfying  the  boundary  conditions.  By  (19),  W/K  is  an 
increasing  function  of  y  for  any  fixed  5,  so  that  a)2  is  also  an  increasing 
function  of  the  adiabatic  index  y.  An  incompressible  medium  (  y  +  °°)  is  thus 
more  stable  than  any  with  finite  y. 

Newcomb has  shown  how  this  result  can  be  extended  to  the  higher  modes 

of  the  system.  If  Z  is  the  subspace  spanned  by  £L,  L  ,  the 

n  wj  ~rl  ^n- 1 

eigenvectors  corresponding  to  w2,  oj2,  ...»  w 2  ^ ,  then  the  energy  principle 


for  the  next  eigenvalue  takes  the  form 


uj2  =  min  (W/K) , 

S  e  c( r  ) 

n 


(22) 


where  C( Z)  is  the  orthogonal  complement  of  Z.  Let 


F  (I)  =  min(W/K),  (23) 

E  C(  S) 

so  that  io2  =  F( Z  ).  The  subspace  Z  is  distinguished  from  all  others  of 
n  n  n 

dimension  n  as  that  which  maximizes  the  function  F  (because  it  is  spanned  by 

the  eigenvectors  with  the  n  lowest  eigenvalues).  Hence 

oj2  -  max  min  (W/K),  (24) 

n  dim(  Z)  «  n  £  e  c(  Z) 

from  which  it  follows  that  w2  is  an  increasing  function  of  y  for  all  n  =  1, 
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We  saw  (in  the  preceding  section)  that  the  degree  of  compressibility 
plays  no  role  in  determining  the  threshold  for  instability.  The  present 
result  implies  that,  of  two  otherwise  identical  unstably  stratified  fluid 
systems,  the  more  compressible  one  has  the  larger  growth  rate.  In  the  sequel 
we  illustrate  this  conclusion  by  calculating  growth  rates  in  some  situations 
where  analytic  solutions  are  possible. 
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Rayleigh-Taylor  Instability  at  the  Interface  between  Two  Isothermal  Layers 
An  ideal  polytropic  fluid  in  a  constant  gravitational  field  evolves  in 


time  according  to  the  system  of  equations 

3p 


8t 


+  V*pv  =  0; 


3v 


p(-^"+  v*7v)  +  7p  +  fig  =  0; 

3p 


(25) 

(26) 

+  v*7p  +  *ypV»v  =  0.  (27) 

at  *** 

The  condition  for  a  stationary  equilibrium  (v  =  0)  is  expressed  by  Eq.  (1). 
Let  us  leave  p  and  p  otherwise  unspecified  for  a  moment  and  suppose  that  this 
state  is  subjected  to  an  infinitesimal  perturbation  defined  by  the  local 
displacement  £(x,y,z,t)  of  an  element  of  fluid.  Using  the  subscript  1  to 
distinguish  perturbed  quantities,  we  have  for  the  velocity 


Z'  It 

so  that  the  perturbed  density  satisfies 


3pj 

IF" 


V,pv  =  -  v(pO, 
1  3t  ^ 


(28) 


(29) 


whence  Pj  is  given  by 


P,  -  -  7-(  pO  (30) 

1 

plus  a  time-independent  quantity  which  can  be  set  equal  to  zero.  Likewise, 

the  perturbed  adiabatic  law 

9p1 

-r— -  +  V  *Vp  +  'TpV'V  =  0  (31) 

d  L  ***1 

has  the  solution 


p,  =  -*Yp7*^  -  ?*7p. 

Substitution  of  (30)  and  (32)  in  the  perturbed  momentum  equation 

p|^+  +  Pjg,-  0 

yields 

;}2  r 

p4rT~  V(^pV*C  +  C-Vp)  -  gV.(pC)  =  0. 


(32) 


(33) 


(34) 
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Assuming  that  the  perturbed  quantities  vary  sinusoidally  with  frequency  w  as 
functions  of  time,  we  can  rewrite  (34)  in  the  form 

-102S  +  (Y-l)go  -  7a  +  7(g-0  =  o,  (35) 

where  a  =  It  is  this  equation  which  must  be  solved,  subject  to  boundary 

conditions,  in  order  to  determine  the  eigenvalues  to2  and  the  spatial 
dependence  of  the  eigenvectors 

Evidently  the  simplest  choice  of  the  basic  state  is  isothermal, 

p/ p  =  c2  =  const,  (36) 

which  leads  to  an  equation  for  £  all  of  whose  coefficients  are  constant.  If 
g  is  directed  downward,  i.e.,  in  the  negative  y  direction,  we  then  have 

p(y)  =  Pq  exp  (-gy/c2).  (37) 

As  will  be  seen,  the  specification  (36)-(37)  for  the  basic  state  results 
in  eigenfunctions  which  are  likewise  exponential  in  y,  and  yields  an 
algebraic  dispersion  relation*  Any  other  choice  gives  rise  to  transcendental 
functions  which  must  be  evaluated  at  ordinary  (nonsingular)  points  when  the 
boundary  conditions  are  imposed,  greatly  complicating  the  form  of  the 
dispersion  relation. 

For  our  present  purposes  it  suffices  to  consider  a  piecewise  isothermal 
state  with  just  two  regions  (Fig.  1).  We  take  the  interface  separating  them 
to  be  at  y  =  0,  i.e.  ,  coinciding  with  the  x-z  plane.  To  distinguish  between 
the  regions  we  will  label  all  quantities  belonging  to  the  lower  one  with 
bars.  For  the  sake  of  generality  we  allow  the  adiabatic  exponent  to  vary, 

y  *  y. 
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the  scale  height  in  the  upper  region  and  the  gravitational 


With  this  choice,  the  coefficients  of  Eq.  (35)  become  constants. 
Assuming  harmonic  dependence  in  the  transverse  plane,  with  the  x  axis  chosen 
parallel  to  the  wave  vector  k,  we  can  look  for  solutions  which  are 
exponentials  in  y: 

£(x,y,t)  =  (e  A  +  e  B)  exp [i(kx-ojt)  -  yy],  (38) 

fvr-  ywX  r~y 

A,  B,  and  y  constant.  Thus  (35)  becomes 

-o)2A  -  ikc2(ikA  -  yB)  +  ikgB  =0,  (39) 

-w2B  +  (y-l)g(ikA  -  yB)  -  yc  2(ikA  -  yB)  -  ygB  =  0.  (40) 

Setting  the  determinant  of  (39)-(40)  equal  to  zero  yields  a  quadratic  for  y. 
The  condition  that  the  solution  be  well  behaved  as  y  +°°  selects  the  root 

y  =  {-yg  +  [  y2g2  -  4  w2yc2  +  4  y2k2c4  -  4y(y-l)g2k2c2/u)2]1/2}(2yc2)~1.  (41) 

Similarly  we  obtain  a  quadratic  for  y  in  the  lower  half-plane;  the  condition 

that  the  eigenfunctions  vanish  at  y  =  -00  yields 

"u  =  {-~Yg  -  ["Y^g2  -  4w2yc2  +  4"y2k2c4  -  4"y(ry-l)g2k2c2/ w2]  1/2}(2~yc  2)“  (42) 

The  kinematic  boundary  condition  at  the  interface  reduces  to  continuity 

of  the  vertical  component  of  i.e., 

A  =  A.  (43) 

The  dynamic  boundary  condition  requires  that  p1  +  be  continuous  at 

y  *  0,  whence  by  Eq.  (32) 
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YP ( 0)  o  -  yp(0)  Q 


(44) 

Since  p(0)  =  p(0),  we  can  combine  (43)  and  (44)  as 

y(ikA  -  uB)/A  =  y  (ikA  -  jjB)/A.  (45) 

Substituting  B/A  from  (39)  or  (40)  and  writing  \i  in  terms  of  u)2  by  means  of 

(41) ,  and  performing  the  analogous  operations  with  the  barred  counterparts  of 
these  equations,  we  find  from  (45)  a  relation  determing  w2.  When  y  =  y  it 
becomes  formally  identical  with  the  wave  dispersion  relation  given  by 
Tolstoy10,  and  for  the  special  case  y  =  y  =  1  it  reduces  to  the  one  given  by 
Blake. 8 

The  treatment  can,  however,  be  carried  a  step  further.  Using  the 
interactive  symbolic  manipulation  system  MACSYMA,  we  can  transform  this 
equation  into  an  algebraic  equation  in  Z  =  u)2/kg.  This  is  done  by  squaring 
the  equation  twice  to  eliminate  the  square  roots  appearing  in  Eqs.  (41)  and 

(42)  and  factoring  the  result.  The  physical  root  is  found24  to  satisfy  the 
quartic 

D'  2  Z4  -  2DDf  SZ 3  +  (D2S 2  +  2DD!  -  Df  2)Z2  -  2  D2(S  -  S'  )Z  -  D4  =  0,  (46) 

where 


D  = 

k(c 2  -  c2)g-1; 

(47) 

S  = 

k(c  2  +  c  2)g- 

(48) 

D!  = 

k(c2/  T  -  c2/  T)g“  1; 

(49) 

S’  = 

k(c 2/  y  -  c2/y)g-1; 

(50) 

Evidently  Eq.  (46)  always  has  one  negative  root,  which  for  D  <  0  is  found 
numerically  to  satisfy  Eq.  (45).  This  solution  can  be  exhibited  by  applying 
the  general  procedure  for  solving  a  quartic,  but  the  result  is  far  too 
cumbersome  to  be  useful.  Instead  we  look  at  some  limits  and  special  cases  of 
physical  interest. 
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First  let  y  =  y,  so  that  both  regions  contain  fluids  with  the  same 
compressibility  properties.  Then  Eq.  (46)  becomes 

Z4  -  2ySZ3  +  (  y2S 2  +  2y  -  1)Z2  -  2  y(  Y  -  1)SZ  -  Y2D2  =  0.  (51) 


In  the  limit  y  +  <»,  the  solution  of  Eq,  (51)  associated  with  the  instability 
satisfies 

S 2Z 2  -  2SZ  -  D2  =  0,  (52) 

whose  negative  root  is  given  by 

Z  =  [1  ±  (1  +  D2)  1/2]S_1.  (53) 


For  negative  values  of  D  the  lower  sign  in  Eq,  (53)  yields  a  solution  of  Eq. 
(45),  as  confirmed  by  numerical  evaluation.  When  we  take  kc2  »  g,  kc2  >>  g 
(wavelength  short  compared  with  both  scale  heights),  we  recover  the  usual 
dispersion  relation  for  the  Rayleigh-Taylor  instability  at  an  interface 
between  two  uniform  incompressible  media,  viz., 


pn  ~  pn 


kg  P0  +  P0 

At  long  wavelengths  (k  +0),  Eq.  (53)  goes  over  to 


(54) 


o2  =  -k 2 


(c2  -  c2)  2 


(55) 


2(c  2  4*  c  2) 

displaying  the  effect  of  the  spatial  dependence  of  the  unperturbed  state.  We 
thus  recover  the  incompressible  result,  as  expected.  Figure  2  illustrates 
the  approach  to  this  limit.  It  shows  plots  of  u)2/kg  as  functions  of  k 
obtained  by  solving  Eq.  (46)  numerically  for  various  choices  of  y  between  1 
and  °°,  assuming  the  unperturbed  states  shown  in  Fig.  1.  As  can  be  seen, 
finite  compressibility  increases  the  growth  rates,  the  relative  effect  being 
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greatest  at  long  wavelengths.  When  k  +  0  then  S,  D,  S!  and  D1  all  become 
small  and  Eq.  (46)  reduces  for  finite  y,  y  to 

D' (2D  -  D')Z2  -  2D2( S  -  S')Z  -  D4  =  0,  (56) 

whence  u)  is  proportional  to  k.  For  y  =  y  Eq.  (56)  yields 

Z  =  Y-1)S  -  [ (  Y-l)  2  +  (2  Y-1)D2]  1/2}.  (57) 

When  y  =  1  this  becomes 

oo2  =  k2(c2  -  c2).  (58) 

This  is  to  be  compared  with  the  corresponding  incompressible  result  given  in 
Eq.  (55).  On  the  other  hand,  for  short-wavelength  perturbations  (k  ■*<»), 

Eq.  (46)  reduces  to 

S2Z2  -  D2  =  0,  (59) 

whose  solution  is  identical  with  Eq.  (54). 

Another  interesting  limit  is  that  in  which  the  density  of  the  upper 
medium  becomes  infinite,  so  that  c  =  0.  One  of  the  extraneous  roots  factors 
out  of  Eq.  (46),  which  then  reduces  to  a  cubic, 

Z3  +  (2y  -  l)z  -  k--Y-  C-2  (Z2  -  1)  =  0.  (60) 

g 

Equation  (60)  holds  even  if  y  +  00  in  such  a  way  that  yc2  (the  square  of  the 
sound  speed)  remains  finite.  If  instead  we  assume  that  the  density  in  the 
lower  region  vanishes  (i.e.,  c2  +  °°) ,  then  the  dispersion  relation  is  even 
simpler,  becoming 

Z  «  -1,  (61) 

which  coincides  with  the  incompressible  result.  This  is  consistent  with  the 
behavior  shown  in  Fig.  2,  which  indicates  that  as  pQ  decreases  (for  fixed 
pQ),  the  difference  between  compressible  and  incompressible  growth  rates 
becomes  less  pronounced. 
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0.3 


Figure  2 

Dimensionless  squared  growth  rate  -co2/kg  vs  wavenumber  k  for  the  two 
basic  states  shown  in  Fig.  1,  with  the  same  choice  of  units.  The  adiabatic 
index  y  in  both  regions  is  taken  to  be  1,  5/3  or  °°,  as  indicated  by  the 
label.  Note  that  the  curves  asymptotically  approach  the  value  (p  -  p)/ 

(p  +  p),  equal  to  0.333  and  0.818,  respectively. 


22 


If  we  assume 


V"  PqY, 

then  D1  vanishes  identically  and  Eq.  (46)  simplifies  to 

S2Z2  -  2(S  -  S'  )Z  -  D2  =  0, 


whence 


Finally,  if 


s  -  S'  -  [(S  -  S')2  +  P2S2]  1/2 

2S-2 


P0  P0  ’ 

so  that  c2  =  c2  and  D  =  0,  then  even  for  "y  *  y 

Z  =  0, 

i.e.,  the  perturbations  are  marginally  stable. 


(62) 

(63) 

(64) 

(65) 

(66) 


23 


Use  of  Piecewise  Isothermal  Profile  to  Approximate  an  Arbitrary  Equilibrium 
More  generally,  we  can  specify  an  equilibrium  state  consisting  of  N-l 
slabs  of  finite  thickness  sandwiched  between  two  semi-infinite  regions,  with 
density  profiles  in  the  various  regions  of  the  form 

pj(y)  =  pj  exp  (-gy/c|),  (67) 

pressure  profiles  given  by 

Pj(y)  =  Pj(y)  c|,  (68) 

and  adiabatic  indices  Yj  ,  j  *  0 ,  1 , . . .  ,  N.  We  take  the  interface  separating 
layer  j  from  layer  j+1  to  be  at  y  =  y^ ,  j  =  0,  1,  . ..,  N-l,  with  y^  «  0.  At 
each  interface  the  changes  discontinuously  but  the  pressure  is  continuous,  so 
that 

pjcj  exp  (-gyi/cj)  =  pj+1cj+1exp  ( -gy.j /cj+1) ,  (69) 

j  «  0,  1,  . ..,  N-l.  Evidently  such  a  piecewise  isothermal  state  can  be  made 
to  approximate  an  arbitrary  ideal  hydrostatic  equilibrium  state  as  closely  as 
desired  if  the  number  of  interfaces  N  is  allowed  to  increase  without  bound. 

It  is  thus  analogous  to  the  piecewise  isopycnic  (constant-density)  model  used 
by  Mikaelian2^  to  approximate  an  arbitrary  incompressible  equilibrium  state. 
Following  the  treatment  employed  in  the  previous  section,  we  seek  a 
solution  for  the  perturbed  displacement  in  the  form 

Jo  -  (^V-yV  exp[i(kx-ut)-y*y],  (70) 

+  + 

Jg  =  I  C^Aj  e  +v%B1~e~V1jy)  exp[i(kx-0Jt)  ] ,  (71) 

j  =  1,  2,  ...  N-l,  and 
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(72) 


4  =  +^yV  -  V1* 

where 

y.”  =  {-Y.g  ±  [Y?g2  -  4cj2y.  c.2  +  4Y?k2c.4 
J  J  J  J  J  J  J 


“  ^Y.(  Y.-Dg^^?/  u2]  1/2}(2y.c  2)~  1  .  (73) 

J  J  J  J  J 

+  + 

Evidently  Eqs.  ( 7 0)  —  ( 72)  introduce  4N  unknown  quantites  A ^  ,  B_.  .  The 
boundary  conditions 


+  -v.y 


l  A .  e  J  J  =  l  A,  . .  e 


±  -^4.1  y. 


j+i'j 


(74) 


+  J 


+  J+1 


and 


± 


YjJ  (lkAj  yjBj)e  j  Yj+ii:  (ikAj+i  Mj+iBj+i)e 


(75) 


j  =  0,  1,  N-l,  provide  2N  linear  relations  among  these*  (Note  that 

—  —  -f  -f  ±  ± 

A^  =  B^  =  A^  =  B^  =  0. )  We  can  eliminate  the  B^  in  favor  of  the  A^  using  the 


analog  of  (39), 


-u>2  Aj  -  ikc^2  (ikA"  -  Uj  B  )  +  ikgB^  =  0, 


(76) 


leaving  2N  linear  homogeneous  equations  in  the  2N  quantities  A^ .  The 


dispersion  relation  giving  u)2  as  a  function  of  k  is  then  obtained  by  equating 
to  zero  the  determinant  of  this  system. 
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The  resulting  equation  must  be  solved  numerically.  It  is  found  that  the 

number  of  roots  associated  with  gravity  modes  equals  N-l,  the  number  of 

interfaces.  For  short  wavelengths  (k  +  *>) ,  one  mode  is  localized  at  each 

interface  position  y  and  is  stable  or  unstable  according  as  or 

pio+1  >  i-  At  longer  wavelengths  the  identity  of  various  modes  becomes 

obscure,  and  the  distinction  between  Rayleigh-Taylor  and  convective 

instablity  may  be  blurred.  Many  interesting  limiting  cases  can  be 

distinguished,  e.g.  ,  y^  +  °°,  j  =  0,  1,  ...,  N  (incompressible  fluid); 

0  N 

Pq  +  00  (solid  lower  boundary);  +  0  (free  upper  surface),  etc.  Of  course, 
it  might  be  argued  that  it  is  just  as  easy  (in  both  incompressible  and 
compressible  cases)  to  approximate  the  differential  equation  for  £  by  finite 
differences  and  obtain  the  spectrum  using  a  standard  eigenvalue-solving 
routine. 
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Uniform  Self -Similar  Implosions  and  Expansions:  Formulation 

Many  of  the  problems  in  which  compressible  fluids  are  subject  to 
Rayleigh-Taylor  instability  involve  nonsteady  basic  states  (e.g. ,  laser 
pellet  implosions,  the  outward  motion  of  gas  following  a  stellar  explosion). 
Usually  in  such  problems  even  the  unperturbed  motion  can  only  be  treated  by 
solving  the  fluid  equations  numerically  using  a  code  with  one  spatial 
variable.  To  analyze  perturbations  with  angular  dependence  requires  a  two- 
or  three-dimensional  code.  A  compromise  approach  involves  linearizing  the 
perturbed  equations,  expanding  them  in  cylindrical  or  spherical  harmonics, 
and  advancing  the  perturbed  radius-dependent  amplitude  functions  in  parallel 
with  the  variables  describing  the  basic  state.17 

However,  there  exists  a  class  of  nontrivial  ideal  compressible  flows 
which  are  sufficiently  symmetric  that  both  the  unperturbed  and  perturbed 
equations  can  be  solved  analytically.  These  are  the  uniform  self-similar 
solutions  studied  by  Sedov. 25  They  can  readily  be  derived  as  follows. 27 
We  rewrite  Eqs.  (25)-(27)  in  the  form 


p  +  p7*v  =  0: 

(77) 

P  +  Vp  =  0; 

(78) 

(pp“Y)  *  =  o, 

(79) 

. 

where  the  raised  dot  (  ) 

denotes  a  total  time  derivative. 

In  a  system  with 

planar,  cylindrical,  or 

spherical  symmetry  (v  =  1,  2,  3, 

respectively) , 

Eqs.  (77)— (78)  become 

p  +  pR  V  -|R  (RVv) , 

(80) 

P  V  + 

(81) 

27 


For  uniform  self-similar  flow  there  is  a  function  f(t)  such  that  for  an 
arbitary  fluid  element  whose  position  at  t  «  0  was  r,  the  position  R  at  time 
t  satisfies 

R  -  rf(t),  (82) 


with  f ( 0)  =  1  and  f ( 0)  =  0.  The  continuity  equation  (80)  then  yields 

p(r,t)  =  pn(r)  f_V, 


(83) 


and  hence  from  the  adiabatic  law  (79) 

p(r,t)  =  p0(r)f_v,=  s(r)P(Jf' 
where  we  have  introduced  the  entropy  function  s(r)  =  p^ 


-Y 


(84) 

If  we  specify 


an  initial  density  profile  p^(r)  on  some  interval  r^  <  r  <  rQ,  then 


substitution  in  Eq.  (8)  results  in  separation  into  a  spatial  part 

dP , 


0 


dr 


-2 


which  determines  the  pressure  Pq,  and  a  time-dependent  part 


f®1-1  f  =  +t"? 


(85) 


(86) 


where  a=  v( y- 1 ) ,  determining  f.  In  Eqs.  (85)-(86)  x  is  a  separation 
constant  with  units  of  time;  the  upper  (lower)  sign  corresponds  to  implosion 
(expansion).  If  the  pressures  on  the  inner  and  outer  surfaces,  p(r^,t)  *= 

—  VY  —  VY 

Po(rf)f  and  p(rQ,t)  =  p^(rQ)f  ,  are  nonvanishing,  they  must  be  balanced 
by  an  equal  pressure  applied  to  the  shell.  It  is  difficult  to  imagine  how 
this  might  be  realized  in  practice,  so  we  will  assume  Pp(r^)  =  0  for 
implosions  and  p^(r^)  =  0  for  expansions. 

A  quadrature  can  be  performed  on  (86),  with  the  result 


2  *2 

xzfz  =  +2  ftif 


(87) 
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when  y  =  1,  and 


T2f2=  t!  (l-f-CX)  (88) 

otherwise.  When  a  =  2  (corresponding  to  y  =  3,  2,  5/3  for  v=  1,  2,  3, 
respectively),  Eq.  (88)  can  be  integrated  directly  to  give 

f(t)  =  (1  +  t2/x2)l/2.  (89) 

For  other  values  of  y,  f  is  most  conveniently  found  by  numerical  means.  In 
the  case  of  implosions,  f  vanishes  at  a  time  t^  given  by 


1/2  1 

VT-(f)  / 


df 


0  (f"a-  1) 


rjr_  1/2  rq/a  +  i/2) 

1/2  "  2a  r(  1/ a  +  1) 


(90) 


which  decreases  monotonically  as  a  function  of  ou 

To  study  the  stability  of  uniform  self-similar  flows  under  small 
perturbations,  we  must  solve  the  linearized  form  of  (78)  for  the  first-order 
displacement 


P(  5  -  S-V_R)  -  R7_  -  V  (pV  +  5-V  p)  «  0. 

K  ^  K  K  ^  K 


(91) 


Note  that  (91)  is  identical  with  (34),  except  that  g  is  replaced  by  -R.  Like 
the  unperturbed  momentum  equation,  (91)  is  separable  in  Lagrangian  variables. 

Substituting 


£(r,t)  «  5(r)T(t), 


(92) 


we  have,  on  writing  7=7, 


(  X  -  1)  5  +  (  Y  ~  1)  r  a  ±  t  (  YPn/  Pn)  Va  +  7(r  •£)  =  0 


(93) 


[cf.  Eq.  (35)],  and 


2  j.cri-2" 

t  f  T  =  +  XT, 


(94) 
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where  X  is  the  new  separation  constant,  obtained  as  an  eigenvalue  in 
connection  with  the  solution  of  Eq.  (93).  Once  X  is  known,  (94)  can  be 

converted  into  a  hypergeometric  equation  in  the  new  variable  x  =  1  -  f  a  and 
solved  to  give 

T(t )  =  T(0)F(a ,  b;  1/2;  x)  +  tT(O)  ( +2x/ a) zF(a  +  1/2, -b  +  1/2;  3/2;  x) 


(95) 

where  F  is  the  hypergeometric  function. 23 

Here 

a  =  (  a  +  2  + 

A) /A  ct, 

(96a) 

b  =  (  a  +  2  - 

A) /A  a, 

(96b) 

n 

with  A  =  [  a  +  2)  -  8aX]1/2.  When  y  1 ,  Eq.  (95)  goes  over  to 


T(t)  =  T(  0)  $(  X/  2 ;  1/2;  tof)  6  tT(0)  (  +2  tof )  1/2  $<  A/2  +  1/2;  3/2;  tof),  (97) 

where  $(a;  b;  x)  is  the  Kummer  function. 29 

Since  T(t)  is  not  exponential,  we  must  decide  what  we  mean  by 
instability.  A  perturbation  is  defined  to  be  unstable  if  the  associated 
time-dependent  factor  satisfies 

tliPco  lT(t)  |/f<t)  -  -»  (98) 

and  stable  if  the  limit  of  this  ratio  is  finite.  This  is  equivalent  to 
saying  that  a  perturbation  is  unstable  if  and  only  if  its  amplitude 
eventually  becomes  infinitely  larger  than  the  radius  of  the  unperturbed 
state. 
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Stability  of  Uniform  Self-Similar  Implosions 

We  begin  by  considering  imploding  systems.  As  f  +  0,  both  terms  in  Eq . 

(95)  approach  asymptotic  forms  containing  terms  proportional  to 

f  (^^-A)/^  with  the  lower  sign  this  expression  diverges  whenever 
a  +  2  <  A,  i.e.,  X  <  0.  Since  the  condition  for  A  to  be  real  is  that  X  be  no 

greater  than  (a  +  2)2/8a,  whose  minimum  value  as  a  function  of  a  is  unity,  we 

see  that  X  <  1  is  always  sufficient  to  make  T/f  diverge.  Thus  X  >  1  is  the 

stability  criterion  for  uniform  self-similar  implosions.  If  A  is  imaginary, 

T/f  still  diverges  when  a  <  2,  i.e.,  y  <  1  +  2/  v.  Elsewhere30  I  have 

presented  a  simple  argument  involving  conservation  of  wave  action  to  show 

that  this  describes  sound  wave  amplification  as  a  consequence  of  the 

geometric  properties  of  the  implosion. 

Taking  the  scalar  product  of  (93)  with  £  and  introducing  notations  for 

the  transpose  V£+  and  the  curl  w  =  Vx£,  we  can  multiply  through  by  Pq  and 

integrate  to  obtain  an  energy  principle:31 

t"2  X/dVpnS2=  JdVpn[(Y-l)  a2  +  VS:V£*--  a)2]  +  y  JdSpnan.?,  (99) 
VUVU  s - 

where  V  is  the  volume  and  S  the  surface  of  the  shell,  and  n  is  the  unit 
vector  defined  so  as  to  point  away  from  the  shell  on  both  inner  and  outer 
surfaces.  The  expression  multiplying  X  and  the  first  two  terms  in  the  volume 
integral  on  the  right-hand  side  are  manifestly  nonnegative.  From  this  we  see 
that  relative  instability  (  X  <  1)  can  only  result  if  u>  *  0  or  if  the  surface 
integral  is  nonvanishing.  The  latter  is  the  case  whenever  a  perturbation 
exists  at  a  point  where  the  density  changes  discontinuously ,  e.g. ,  at  r  =  r  , 
r  =  r^,  or  an  internal  density  jump.  The  external  pressure,  which  enters 
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the  model  as  a  boundary  condition,  produces  an  inward  acceleration.  There  is 
thus  an  effective  gravity  in  the  outward  direction.  Hence  we  anticipate  that 
a  Rayleigh-Taylor  instability  should  occur  localized  at  the  outer  surface. 
This  case  was  first  studied  by  Kidder3^,  who  assumed  y  =  5/3.  It  turns  out 
we  can  solve  for  £  provided  the  perturbation  wavelengths  also  expand  or 
contract  self -similarly  in  time,  i.e.,  whenever  they  do  not  introduce  a 
definite  length  scale  into  the  problem.  This  means  that  in  spherical 
geometry  ( v  =  3)  there  is  no  restriction  on  the  form  of  the  perturbation;  in 
cylindrical  geometry  (v  =  2),  however,  we  must  have  =  0;  and  no  general 
solution  is  possible  in  planar  geometry  (v=  1). 

Operating  on  (93)  with  the  divergence  and  with  the  curl,  we  get  two 
equations : 


[  \  +  v(  y-l)  +  1  ]o  ±  x  V*( -  Va)  +  yr  »Va  =  r  • 

po  *~  ~~ 


Vxijj 


(100) 


and 


(  X-l)  w  +  r  xVa  ±  :2ip  7(  p  *)  x  Va  =  0.  (101) 

u  u 

Let  us  suppose  that  a> =  0.  It  follows  from  ( 1 00)— ( 101)  that  a  also  vanishes. 
If  that  happens, 


5  *=  7i|>,  (102) 

where  the  potential  ip  satisfies  Laplace’s  equation 

V2tJj  =  0.  (103) 


The  general  solution  of  ( 103)  is 

<K  r  ,<())“(  t+r  +  1^_r  )  e 


(104) 


in  cylindrical  coordinates  (assuming  -r—  =  0),  and 

oZ 

<K r,  9,  4>)  =  (  r  %  +  r-"'-1)  Y ^(  0,  <j>)  (  105) 

in  spherical  coordinates.  The  ^ ,  are  constants,  and  the  Y «  (  0,  4>)  are 

x  Jan 
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spherical  harmonics.  Evidently  the  first  term  in  (104)  and  (105)  corresponds 
to  a  solution  localized  at  the  outer  surface  of  the  shell,  and  the  second  to 
one  localized  at  the  inner  surface,  so  we  set  «  0.  Substitution  in 
Eq.  (93)  then  yields 

(X-l)5+  V(r  •£)  =  v[(X-l)  <|>+  r«V+]  -  V[  (  X-l+Z)  * r  l)  =  0,  (106) 

****■  *****  *****  T 

whence 


X  -  -A  +  1.  (107) 

For  £=0,  X-  1  and  Eq.  (94)  reduces  to  Eq.  (86),  showing  that  the 
perturbations  are  marginally  stable.  For  all  l  >  0,  the  limiting  form  of  T/f 
diverges  as 

T  _  f  {®-2-[(»-2)*f8al]l/2}  (10g) 

when  f  +  0.  Evidently  the  magnitude  of  the  exponent  in  (108)  increases  with 
both  Z  and  a=  v(y-l)*  Thus,  in  contrast  with  the  case  of  static  equilibria 
considered  previously,  compressibility  appears  to  be  somewhat  stabilizing. 
However,  it  must  be  noted  from  Eq.  (86)  that  as  a  increases,  the  motion 
becomes  increasingly  stiff  and  so  the  effective  gravity  also  increases, 
rendering  comparisons  difficult. 

The  problem  of  perturbations  with  u>  *  0  is  treated  elsewhere30-31;  it  is 
completely  analogous  to  that  in  the  case  of  expanding  solutions  to  be 
discussed  shortly.  Book  and  Bernstein33  and  Han  and  Suydam3l+  have  treated 
the  stability  of  imploding  cylindrical  systems  in  detail. 
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Stability  of  Uniform  Self-Similar  Expansions 

Here  we  consider  expanding  systems,  in  which  f  +  «>.  (Sedov 
distinguishes  a  third  class  of  trajectories  in  which  f  varies  between  0  and  00 
with  no  turning  point;  we  do  not  treat  them  here.) 

Bernstein  and  Book35  and  Han35,  who  found  the  general  solutions  for  the 
time  dependence  of  arbitrary  perturbations  in  spherical  and  in  cylindrical 
geometry,  respectively,  both  assumed  that  the  unperturbed  states  were 


—  v 

homentropic  (p p  independent  of  radius).  In  both  geometries  the  only 
instability  was  a  Rayleigh-Taylor  mode  localized  at  the  inner  surface, 

where  the  driving  pressure  acts.  Here  we  analyze  a  class  of  states, 

parametrized  by  the  adiabatic  index  and  a  shape  parameter,  which  relax  the 

requirement  of  uniform  entropy  distribution. 3 7 

Suppose  the  initial  density  profile  is  given  by: 

Pq(t)  =  p(l-r2/r2)K,  (109) 

y  and  p  constant.  Then  from  Eq.  (85), 

pQ  (r)  =  p(l-r2/r2)  *fl,  (110) 

where 

p  =  1r2/ 2(  k+1)  t2,  (111) 

and 

s(r)  =  PqPq  =  (P/P  )(l-r  /rQ)  (112) 

The  pressure  vanishes  at  the  outer  radius  of  the  shell,  as  does  the  density 
provided  k  >  0.  At  the  inner  radius  the  imposed  (driving)  pressure  must  have 
the  form 


✓V  i  2  ,  2v  kH-1  --  vy 
P.  -  P(l-^/r0)  f  . 

unless  jy  =  0.  In  this  model  the  temperature  0  =  p/ p  always  decreases 


(113) 
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quadratically  as  a  function  of  radius: 

9(r,t)  =  0(l-r2/r^)f"v(T"1),  (114) 

where  9  =  r2/2(KH-l)i? 

As  in  the  implosion  case,  there  is  an  incompressible  irrotational 
perturbation  mode  satisfying  Eqs.  ( 102) — ( 103) •  The  solution  (which  is 
independent  of  the  shapes  of  the  density  and  pressure  profiles)  is  given 
again  by  (104)  or  (105).  This  time,  though,  it  is  the  mode  corresponding  to 
which  is  unstable,  giving  rise  to  eigenvalues  X  =  i  +  1  (cylindrical 
geometry)  or  X  =  1+2  (spherical  geometry).  Using  standard  formulas 28  to 
evaluate  hypergeometric  functions  of  unit  argument,  we  find  from  (95)  that  as 
f  -►  », 

T  + _ r(  1/2)  rpToQ  T(0)  +  (2/ q)  1/2  1X3/2)  r(l/a)iT(0)  .(n5) 

f  r(l  +  (2+A)/4a)  r[i  +  (2-A)/4a]  r[|  +  (2+ A) /A  a]  T[j  +  (2-A)/4  a] 

Although  this  limit  is  finite,  for  large  values  of  X  (large  £)  the  constants 
are  found  from  Stirling1 s  formula  to  grow  exponentially: 

Z  ~  T(  1/  a)  exp  [ tt(  X/2  a)  ]_ 1 

f  '  2»1/2  (V2«) 

Figure  3  displays  the  late-time  asymptotic  behavior  of  these  solutions  as  a 
function  of  the  eigenvalue  X  for  a=  1,  2,  3,  and  °°. 

To  study  the  modes  for  which  ^co  *  0,  we  use  Eqs.  (100)-(101)  with  and 
Pq  given  by  (109)-(110).  Using  the  lower  sign,  specializing  to  v  *  3,  and 
using  r^  to  scale  jr,  we  get 

(X  +  3y  -  2)o  -  2(-^hT  fV’(1  "  r2)  Vo^l  +  V70  "  (117) 

and 

(X  -  1)M  -  (  Y*  -  1)  rxVo.  (118) 

***  K  +  1 


T(0)  +  ru/°>  exp[«(V2oO>'2]  T  i(0). 
2(  2  ira)  1/2(  X/2  a)  (  2+a)/4  a 


(116) 
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co— i 


0)J/0!*‘X)£  wji 
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(The  corresponding  limit  for  £)  is  infinite. 


Eliminating  between  (117)  and  (118)  and  assuming  separation  of  the  angular 
and  radial  dependence  of  a  by  setting 

a(£)  =  a(r)  Y^(6,  <j>),  (119) 

we  obtain  a  second-order  equation  for  the  radial  factor  o(jr) , 


1  rl-r2  d  ,  ?d  a 
2(  k+1  )  L  v2  dr  ^  dr; 


£(£+l)a  -  2r 


do  r(  <+l- <Y)  £(  £+1) 

dr  L  (  y-i)(  k+1) 


3Y  -  A  +  2  ]  a  =  0. 


(120) 


Rewriting  this  by  means  of  the  substitution  o=  rV  an<l  x  *  r 2,  we  obtain 
the  hypergeometric  equation 

x  ( 1  -  x)  y1  1  +  [c  -  (a+b+1)  x]  y'  -  aby  =  0,  (121) 

where 

£}  =  \  (<  +  £  +  |  ±  [(  k  +  £  +  |)2  -  AK]1/2),  (122a, b) 

c  =  £  +  3/2.  (122c) 


Here 


K 


(<+2)£  K+lr,  ^  «  £(£+l)(ic+l-  <y) , 

2  +  IT”  1 A  "  2  +  3Y  (A  -  1x77  1)  J 


(123) 


The  general  solution  of  (121)  is  28 

y  =  C^F(a ,  b;  c;  x)  +  C^x*  C  F(a  -  c  +  1,  b  -  c  +  1;  2  -  c;  x) ,  (124) 

where  C,,  C0  are  constants.  If  r^  +  0,  the  shell  goes  over  to  a  gas  sphere 
expanding  under  its  own  pressure.  In  this  case  only  the  first  term  in  Eq. 
(124)  is  finite  at  the  origin  and  we  must  set  C2  -  0. 

The  boundary  conditions  are  found  from  the  requirement  that  the 
perturbed  pressure  vanish  at  both  surfaces  of  the  shell.  At  the  inner  suface 
this  implies  a  =  0  or 


37 


(125) 


CjFCa,  b;  c;  r2/r2)  +  C^r^r^1  C 

*F(a  —  c+1,  b  —  c+1;  2-c;  r2/ r2)  =  o. 

Since  the  unperturbed  pressure  vanishes  at  the  outer  surface,  it  is  only 
necessary  that  y  be  finite  at  x  =  1.  The  linear  connection  formulas28  of 
both  terms  of  (124)  contain  a  term  that  diverges  as  (1  -  x)  ^  unless  a  or 

b  is  a  nonpositive  integer.  Thus  we  must  have 

-i  {<  +  l  +  5/2  -  [( k  +  %  +  5/2) 2  -  4  K]  1/2}  =  -n,  (126) 

n  =  0,  1,  2,...  From  (123)  it  follows  that  X  decreases  with  n.  Hence  the 
fastest  growth  (largest  positive  X)  corresponds  to  n  =  0,  which  implies  K  * 

0.  Solving  for  X,  we  finally  obtain  the  dispersion  relation 

,  ,  _  y£(  <  +  2)  +  (  <  +  l)(3y  -  1) 

X  ■  1 - 2(“Tj 

(127) 

4.  {fy£(  <  +  2)  +  (  K  +  l)(3y  -  1)  l2  +  4£(  l  +  1)(  k  +  1  -  <y)  } 1/2 

2(  <  +  1) 

For  the  upper  branch,  X  >  1  for  all  i  >  0,  provided  that 

<  <  l/(y  -  l).  (128) 

From  (112),  we  see  that  this  is  just  the  condition  that  s*(r)  <  0,  i.e.,  that 
s(r)  decrease  in  the  "upward”  direction,  that  directed  opposite  to  the 
effective  gravitational  acceleration  -R.  This  is  identical  with  the 
Schwarzchild  criterion  (18)  derived  for  the  convective  instability  in  static 
media.  Indeed,  an  interchange  argument  along  the  lines  of  that  used  to 
obtain  (18)  has  been  employed37  to  show  that  such  a  uniformly  expanding  fluid 
system  should  be  convectively  unstable  whenever  sf(r)  <  0  holds. 

For  the  y  =  1  case  we  can  redo  all  the  analysis  in  terms  of  confluent 
hypergeometric  functions28  instead  of  hypergeometric  functions,  or  we  can 
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reach  the  same  results  by  formally  letting  y  -►  1  in  the  equations  above.  The 
solutions  are  qualitatively  similar  to  those  for  y  >  1  except  that  now  (128) 
always  holds. 

Motivated  by  experiments  investigating  the  use  of  imploding  cylindrical 
liquid  metal  liners  to  compress  and  heat  plasma38,  Book  and  Bernstein33 
studied  the  Rayleigh-Taylor  instability  on  the  inner  surface  of  a  liner 
during  both  the  implosion  and  rebound  phases,  assuming  y=  1.  Since  both 
terms  in  Eq.  (97)  diverge  the  same  way  at  large  t ,  it  is  useful  to  introduce 
in  their  place  two  new  solutions  which  have  asymptotically  like  ¥,  the 
standard  confluent  hypergeometric  function  of  the  second  kind:2^ 


p*c 


*>) 

=  *K 


l+JO/2;  1/2;  Snf]  +  xf  $(  1+ 1/2 ;  3/2;  taf).  (129a,  b) 


Q*(t)  j 

For  large  arguments  (t  +  +») ,  we  have 


r[(l+£)/2) 


P  V)  ~  (  taf)  (  Z+1)/2 

Q*(t)  ~f(«nf)V2. 


(130a) 

(130b) 


Thus  defined,  0  (t)  has  the  property  of  increasing  monotonically  for  all  t, 

% 

-oo  <  t  <  °°;  at  turnaround  (the  instant  t  =  0  when  f  =  1),  Q  (0)  -  1  also. 

The  only  perturbations  which  are  unstable  both  before  and  after  turnaround 
are  those  whose  time  dependence  is  proportional  to  0  (t).  Figure  4  compares 


the  behavior  of  Q  (t)  for  Z  =  1  and  Z  =  10  with  that  of  f. 
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105 


Figure  4 

% 

Q  (t)  [defined  by  Eq.  (129b)]  for  i  =  1  and  i  =  10,  obtained  by  numerical 
solution  of  Eq.  (94)  for  a=  0  and  X=  i  +  1,  with  f  plotted  for  comparison. 
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Conclusions 


In  the  theory  of  compressible  fluids,  the  Rayleigh-Taylor  instability  at 
a  density  jump,  the  Rayleigh-Taylor  instability  in  a  continuously  stratified 
medium,  and  the  convective  instability  are  close  relatives.  All  are 
gravitational  interchange  modes.  One  can  easily  generate  a  series  of 
examples  which  display  a  continuous  transition  from  one  mode  to  the  next. 

The  energy  principle  applied  to  polytropic  media  shows  that,  by  itself, 
compressiblity  increases  instability.  For  only  a  handful  of  specific 
compressible  states,  however,  is  it  possible  to  actually  calculate  growth 
rates  in  closed  form.  The  only  tractable  equilibrium  states  involve 
contiguous  isothermal  layers  of  fluid  satisfying  the  adiabatic  law  with 
constant  y.  In  the  limit  where  the  density  of  the  lower  layer  vanishes,  the 
growth  rate  reduces  to  the  classical  result  found  for  incompressible  fluids. 

Closely  related  is  the  problem  of  the  stability  of  uniformly  imploding 
or  expanding  shells  driven  by  pressure  applied  at  a  vacuum-material  boundary. 
The  unstable  modes  are  incompressible  and  (if  one  allows  for  the 
nonexponential  time  dependence)  the  growth  rates  are  given  by  the  same 
classical  incompressible  fluid  formula  as  in  the  static  case. 

Because  the  unstable  eigenmodes  are  localized  near  a  density  jump  within 
distances  of  order  k  * ,  one  expects  the  growth  rates  not  to  change  very  much 
when  the  basic  state  is  not  isothermal,  particularly  at  short  wavelengths. 
Thus  dispersion  relation  (46)  is  at  least  qualitatively  correct  most  of  the 
time  in  static  situations. 
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To  determine  growth  rates  with  precision  for  any  but  the  simplest 
piecewise  isothermal  unperturbed  states  or  to  study  the  effects  of  thermal 
conduction,  flow  through  the  interface,  nonlinearity,  etc.,  one  must  resort 
to  computational  means.  Nevertheless,  numerical  solutions  present  their  own 
difficulties.  Validation  against  nontrivial  analytical  solutions  such  as 
those  discussed  in  this  review  is  indispensable  in  the  development  of  any 
code. 
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